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Abstract 

Background: An important question is whether evolution favors properties such as mutational robustness 
or evolvability that do not directly benefit any individual, but can influence the course of future evolution. 
Functionally similar proteins can differ substantially in their robustness to mutations and capacity to 
evolve new functions, but it has remained unclear whether any of these differences might be due to 
evolutionary selection for these properties. 

Results: Here we use laboratory experiments to demonstrate that evolution favors protein mutational 
robustness if the evolving population is sufficiently large. We neutrally evolve cytochrome P450 proteins 
under identical selection pressures and mutation rates in populations of different sizes, and show that 
proteins from the larger and thus more polymorphic population tend towards higher mutational robust- 
ness. Proteins from the larger population also evolve greater stability, a biophysical property that is 
known to enhance both mutational robustness and evolvability. The excess mutational robustness and 
stability is well described by existing mathematical theories, and can be quantitatively related to the 
way that the proteins occupy their neutral network. 

Conclusions: Our work is the first experimental demonstration of the general tendency of evolution to 
favor mutational robustness and protein stability in highly polymorphic populations. We suggest that 
this phenomenon may contribute to the mutational robustness and evolvability of viruses and bacteria 
that exist in large populations. 
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Background 

Proteins are quite tolerant of mutations, al- 
lowing evolution to produce highly diverged 
sequences that fold to similar structures and 
perform conserved biochemical functions [1,2]. 
However, proteins with nearly identical struc- 
tures and functions may differ in their robust- 
ness to mutation [3-5], as well as in their ca- 
pacity to acquire new functions [5]. The fact 
that mutational robustness and evolvability can 
vary among the functionally equivalent proteins 
produced by natural sequence divergence makes 
these properties important hidden dimensions 
in evolution — direct selection for protein func- 
tion is blind to them, yet they can play a cru- 
cial role in enabling future evolution. Whether 
the evolutionary process somehow promotes the 
acquisition of mutational robustness and evolv- 
ability therefore remains a major question [6-8]. 

Previous experiments have identified several 
specific evolutionary conditions that can affect 
mutational robustness. For example, genetic 
complementation decreases the mutational ro- 
bustness of viruses [9] , while high mutation rates 
favor mutational robustness in simulated digi- 
tal organisms [10]. However, theory [11] makes 
the much broader — and heretofore experimen- 
tally untested — prediction that extra muta- 
tional robustness will arise quite generally in 
sufficiently large populations. This prediction 
cannot be understood in the standard frame- 
work of Kimura's neutral theory [12], since one 
of the usual assumptions of the neutral theory 
is that mutational robustness is constant. (Al- 
though Takahata [13] treated the consequences 
of stochastically fluctuating neutrality on the 
molecular clock, he did not describe how muta- 
tional robustness might change systematically 
during evolution.) However, changes in mu- 
tational robustness can be described by envi- 
sioning evolution as occurring on neutral net- 
works, or sets of functionally equivalent pro- 
teins that are connected by single mutational 



steps [14-17]. In a seminal theoretical analysis 
of evolution on neutral networks, van Nimwe- 
gen and coworkers [11] predicted that the extent 
of mutational robustness should depend on the 
degree of population polymorphism. Here we 
briefly summarize their reasoning, since it moti- 
vates our experimental work. We also refer the 
reader to chapter 16 of [8], which contains an 
excellent explanation of the densely mathemat- 
ical work of van Nimwegen and coworkers [11]. 

If an evolving population is mostly 
monomorphic, then each mutation is either 
lost or goes to fixation before another muta- 
tion occurs. The population is therefore usually 
clustered at a single genotype and rarely expe- 
riences mutations, meaning that selection does 
not distinguish between genotypes of different 
mutational robustness. All nodes of the neutral 
network are thus equivalent and will be occupied 
by the population with equal probability [11]. 
On the other hand, a highly polymorphic pop- 
ulation is always spread across many nodes of 
the neutral network. When mutations occur, 
the members of the population at highly con- 
nected nodes have a better chance of surviving, 
causing them to be favored by evolution and 
increasing the average mutational robustness 
[11,17-20]. Specifically, a highly polymorphic 
population occupies each node with a proba- 
bility proportional to its eigenvector central- 
ity [11, 17], a measure of how connected it is 
to other connected nodes (a variant of eigen- 
vector centrality is used by Google's PageRank 
algorithm to rank a webpage's importance in 
the network of internet links [21]). Figure [T]A. 
illustrates how mostly monomorphic and highly 
polymorphic populations are predicted to oc- 
cupy a neutral network. For proteins, changes 
in neutral network occupancy should be man- 
ifested by changes in thermodynamic stabil- 
ity [22], with proteins from highly polymorphic 
populations predicted to be more stable than 
their counterparts from mostly monomorphic 
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populations (Figure [TJ3). Note that the extent 
of polymorphism depends on the product of 
the mutation rate and population size, meaning 
that protein populations of different sizes are 
predicted to evolve to different levels of mu- 
tational robustness and stability even if they 
experience the same mutation rate. 



Results and Discussion 

Design of neutral evolution experiment 

To test whether high population polymorphism 
drives an increase in mutational robustness and 
protein stability, we performed laboratory evo- 
lution experiments on cytochrome P450 pro- 
teins. The basic idea was to neutrally evolve 
P450s under a constant selection pressure in 
populations that were either monomorphic or 
highly polymorphic, and observe whether the 
proteins evolved to different levels of mutational 
robustness and stability. The evolution experi- 
ments started with a P450 BM3 heme domain 
that had been engineered to hydroxylate 12-p- 
nitrophenoxydodecanoic acid (12-pNCA) [23]. 
We imposed the selection criterion that Es- 
cherichia coli cells expressing the P450 had to 
yield lysate with enough active enzyme to hy- 
droxylate a specified amount of 12-pNCA in 40 
minutes. This criterion roughly corresponds to 
the case in which an enzyme must catalyze a 
biochemically relevant reaction at some mini- 
mal level in order for its host to survive. Note 
that other properties such as stability and ex- 
pression level can vary freely, provided that the 
criterion for total activity is met. 

The properties of a neutrally evolving pro- 
tein eventually "equilibrate," much as the prop- 
erties of an isolated physical system under some 
macroscopic constraint tend towards the values 
that maximize the system's internal entropy. 
For proteins, this usually means that stability. 



expression, and activity drift towards their low- 
est tolerable values, since the vast majority of 
random sequences do not encode stable, well- 
expressed enzymes (that is, natural selection 
must work against sequence entropy to maintain 
a functional protein) [22,24]. The initial P450 
had been engineered for maximal activity [23], 
meaning that it was not equilibrated to the more 
mild selection criterion of the experiments. We 
therefore neutrally evolved this initial P450 for 
16 generations, introducing random mutations 
with error-prone PGR and retaining all mu- 
tants that met the selection criterion for total 
activity on 12-pNCA. The procedure used for 
this equilibration evolution was similar to that 
for the polymorphic neutral evolution described 
below. As expected, expression, stability, and 
activity all dropped during the equilibration 
evolution. At the end of the equilibration evo- 
lution, we chose a single sequence as the parent 
for the neutral evolution experiments. The gene 
encoding this parent sequence contained 29 nu- 
cleotide mutations and 13 amino acid mutations 
relative to the initial P450 (Additional Filed]). 

We used this parent gene to begin three 
parallel sets of neutral evolution experiments, 
which we named "monomorphic," "polymor- 
phic," and "unselected" (Figure [2]). The 
monomorphic experiments capture the case 
where the population moves as a single entity, 
the polymorphic experiment captures the case 
where the population spreads across many se- 
quences, and the unselected experiments show 
how the gene evolves in the absence of selection 
for protein function. In all experiments, at each 
generation we used error-prone PGR to intro- 
duce an average of 1.4 nucleotide mutations per 
P450 gene (Table [T|). The mutant genes were 
ligated into a plasmid and transformed into E. 
coli [25], and transformants were selected using 
the plasmid's antibiotic resistance marker. For 
the unselected case, we randomly picked one of 
the mutants, recovered the mutant gene with 
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a plasmid mini-prep, and used this mutant as 
the template for the next generation of error- 
prone PGR. We performed four independent 
rephcates of unselected evolution, evolving each 
for 12 generations. 



For the monomorphic and polymorphic pop- 
ulations, we imposed the selection criterion that 
the P450s hydroxy late 12-pNCA with at least 
75% of the total activity of the original par- 
ent gene. We expressed the P450s in E. coli, 
and then assayed the cell lysates for activity in 
a high-throughput 96- well plate format. The 
total amount of product produced by 80 ^ul of 
clarified lysate in 40 minutes was compared to 
the median of four control wells containing the 
original parent P450 to determine if the mutant 
met the selection criterion. The only differ- 
ence between the monomorphic and polymor- 
phic experiments was the size of the evolving 
populations. In the monomorphic limit, each 
mutation is either lost or goes to fixation before 
the next occurs. We enforced this evolutionary 
dynamic by holding the population size to a 
single protein sequence, similar to the "blind 
ant" random walk of [11]. At each generation, 
we assayed a single mutant. If this mutant met 
the selection criterion, then it was carried over 
to the next generation, corresponding to a neu- 
tral mutation going to fixation. If the mutant 
failed the selection criterion, then the popu- 
lation stayed at the previous sequence for the 
next generation, corresponding to a mutation 
lost to selection. If all of the mutants assayed 
had zero or one mutations, then this proto- 
col would correspond exactly to the equations 
of [11,22]. However, in order to achieve appre- 
ciable sequence evolution on a laboratory time 
scale, we used a mutation rate that sometimes 
produced multiple mutations in a generation. 
We mathematically describe this situation in 
the Mathematical Appendix; here we simply 
note that it is possible to think of each gener- 
ation as introducing a single mutational event 



rather than a single mutation. We performed 
22 independent replicates of monomorphic evo- 
lution, evolving each for 25 generations. 

In the polymorphic limit, the population 
spreads across many sequences. To implement 
this experimentally, we assayed 435 mutants at 
each generation. The selection criterion was 
used to classify each mutant as functional or 
nonfunctional. In neutral evolution, all func- 
tional mutants reproduce with equal proba- 
bility. We therefore pooled equal volumes of 
stationary-phase cultures of each functional mu- 
tant and recovered the pooled genes with a mini- 
prep. The polymorphic evolution experiment 
therefore approaches the equations of [11,22], 
again with the exception that a sequence may 
undergo multiple mutations at a single gener- 
ation. We give the equations describing this 
situation in the Mathematical Appendix. Since 
the population evolves deterministically in the 
polymorphic limit [11, 22], a single replicate 
was performed. Because mutations accumu- 
late more rapidly in the polymorphic experi- 
ments than the monomorphic ones, we evolved 
the polymorphic population for 15 generations 
rather than 25. 



Mutations and mutational robustness 

Figure [3] shows how mutations accumulated dur- 
ing the course of the neutral evolution experi- 
ments (full data are in Table [2] and Additional 
File [2]). Since the unselected protein popula- 
tions evolve without constraint, mutations ac- 
cumulate at the same rate at which they are 
introduced by error-prone PGR, 1.4 nucleotide 
mutations per generation. Because selection 
eliminates mutations that disrupt P450 activ- 
ity, mutations accumulate more slowly in the 
monomorphic and polymorphic populations. 
Mutations accumulate more rapidly in the poly- 
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morphic population than in the monomorphic 
populations. This difference in rates is pre- 
dicted by the equations in the Mathematical 
Appendix to be a consequence of the fact that 
the polymorphic population is more mutation- 
ally robust, and so can tolerate more of the 
possible mutations. 

To test directly whether the polymorphic 
population evolves higher average mutational 
robustness, we measured the fraction of 435 
random mutants that met the selection crite- 
rion. Figure m shows that the polymorphic pop- 
ulation neutrally evolved to a markedly higher 
mutational robustness than the monomorphic 
populations, with 50 ± 2% of the final polymor- 
phic population mutants continuing to function 
versus 39 ± 2% for the final monomorphic pop- 
ulations (Chi-square P-value of 10^^ that these 
values are significantly different). The only dif- 
ference between the two types of populations 
was their size, so evolution has clearly favored 
mutational robustness in the larger and thus 
more polymorphic population. This finding 
represents the first experimental support for 
the prediction that highly polymorphic popu- 
lations evolve excess mutational robustness [11]. 

Theory predicts that the excess mutational 
robustness of a highly polymorphic protein 
population comes from increased protein sta- 
bility [22]. Because the P450 variants unfold 
irreversibly, an equilibrium thermodynamic sta- 
bility AGf cannot be measured. We therefore 
determined stability to irreversible thermal and 
chemical denaturation, two highly correlated 
measures of P450 stability that have previously 
been shown to contribute to mutational ro- 
bustness [5] (see Additional Files [3l [U and [5]). 
Figure [5] shows that proteins from the polymor- 
phic population were in fact more stable than 
their counterparts from the monomorphic pop- 
ulation. We also observed that proteins in the 



polymorphic population tended to accumulate 
to higher levels in E. coli (Figure [5]). Elevated 
expression could be a byproduct of increased 
stability, or it could independently increase mu- 
tational robustness by allowing the proteins to 
better tolerate mutations that decrease codon 
adaptation or reduce folding efficiency. It is pos- 
sible that additional unrecognized biophysical 
factors also contributed to the excess mutational 
robustness of the polymorphic population, but 
no such factors were immediately obvious. 



Interpretation in terms of the P450 neutral 
network 

The higher mutational robustness of the poly- 
morphic population is due to the fact that it 
occupies the P450 gene neutral network dif- 
ferently than the monomorphic populations. 
Measurements from the evolution experiments 
can therefore be used to infer basic proper- 
ties of the underlying neutral network of P450 
genes, as originally noted by van Nimwegen and 
coworkers [11]. In the Mathematical Appendix, 
we derive approximations for the normalized 
principal eigenvalue (z^)oo and the normalized 
average connectivity (i^)o of the neutral net- 
work, where in both cases the normalization 
is obtained by dividing by the network coordi- 
nation number. We obtain (i^)oo = 0.51 and 
{i')o = 0.35 for the P450 gene neutral net- 
work. Our ability to consistently estimate these 
two parameters from four different experimental 
measurements supports the idea that the the- 
ory that we elaborate in the Mathematical Ap- 
pendix appropriately describes the experiments. 
The difference between (z^)oo and {iy)o is a mea- 
sure of the extent to which some P450 neutral 
network nodes have more connections than oth- 
ers. We note that (z^)oo is approximately equal 
to the exponential decline parameter for the 
asymptotic decline in the fraction of functional 
mutants with increasing numbers of random nu- 
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cleotide mutations [3, 26, 27] (see Mathematical 
Appendix). Previous studies looking at this ex- 
ponential decline have reported (i^)oo = 0.7 for 
subtilisin [26], (i^)oo = 0.7 for 3-methyladenine 
DNA glycosylase [27], and {v)^ = 0.7 - 0.8 
for TEMl /3-lactamase [3]. These comparisons 
suggest that P450 has a sparser neutral net- 
work (smaller {v)oo) than these other proteins. 
We suspect, however, that these earlier studies 
(one of which is our own) overestimate (1^)00 
due to insufficient equilibration of the starting 
sequence. We believe that the approach of the 
current work is more accurate for determining 
(i/)(x) because the measurements are made after 
many mutations have equilibrated the initial se- 
quence. This approach could be used in future 
experiments to compare the neutral network 
connectivities of proteins from different fami- 
lies. 



Conclusions 

We have demonstrated that neutral evolution 
favors more mutationally robust proteins when 
the evolving population is highly polymorphic. 
Strikingly, the excess mutational robustness is 
due only to population polymorphism, and so 
will arise in any population of sufficiently large 
size. Our work is the first experimental demon- 
stration of this phenomenon, which is predicted 
to occur quite generally in neutrally evolving 
proteins and nucleic acids [11]. Furthermore, 
we were able to identify one of the biophysical 
factors underlying the increase in mutational 
robustness by showing that proteins from the 
highly polymorphic population are more sta- 
ble. We recognize that evolution in a biological 
context will be more complex. In our exper- 
iments, fitness was the P450's ability to be 
expressed in active form by bacteria grown to 
saturation in an environment with plentiful nu- 
trients. Biological fitness, however, depends on 
numerous additional and subtle effects such as 



the metabolic costs of synthesis or the burdens 
imposed by misfolded molecules. Some muta- 
tions that are neutral in the experiments may 
therefore have deleterious effects in a biologi- 
cal setting [28]. The experiments nonetheless 
capture the overriding constraint that proteins 
retain their biochemical functions. Our success 
in quantitatively explaining the results supports 
the notion that important aspects of protein 
evolution can be described simply in terms of 
mutational effects on stability [22,28]. 

An obvious question is whether evolution 
in nature favors mutational robustness by the 
process we have demonstrated. Whether natu- 
ral populations will neutrally evolve mutational 
robustness depends on whether they are suffi- 
ciently polymorphic, which will be the case if 
the product of their effective population size A'^ 
and per protein per generation mutation rate /x 
is much greater than one [11, 12]. Accurately 
estimating A^^u, which is closely related to the 
widely used parameter 9 in population genet- 
ics, for natural populations is difficult [29, 30] 
(note that since mutational robustness is a 
protein-wide property, the relevant mutation 
rate is per protein, which 10^ to 10^ lar ger 
than the per codon mutation rate). For hu- 
mans and other multicellular organisms, Njj, 
is probably too small [31] for their proteins to 
neutrally evolve mutational robustness. But 
estimates [31, 32] place N jjL k, 10 to 100 for 
typical-length proteins in bacteria, and it is 
probably much higher for many viruses [33,34]. 
It is therefore likely that many viral and some 
bacterial proteins have neutrally evolved extra 
mutational robustness. 



The neutral evolution of protein mutational 
robustness may also contribute to adaptive evo- 
lution. Experiments have shown that extra sta- 
bility increases a protein's evolvability by allow- 
ing it to tolerate a wider range of functionally 



6 



beneficial but destabilizing mutations [5]. A 
similar phenomenon seems to occur in natural 
evolution, where functionally neutral but stabi- 
lizing mutations can play a key role in adaptive 
evolution by counterbalancing the destabilizing 
effects of other functionally beneficial muta- 
tions [35]. Viruses and perhaps bacteria may 
thus benefit from large population sizes and 
high mutation rates that drive an increase in 
the mutational robustness and stability of their 
proteins, which in turn enhances the capacity of 
these proteins to rapidly change their sequences 
and evolve new functions. 



Methods 

Equilibration evolution of the P450 protein 

We began with a 21B3 P450 peroxygenase that 
had been engineered for highly efficient hydrox- 
ylation of 12-pNCA [23] (sequence shown in Ad- 
ditional File [6]) . This P450 was not well equi- 
librated to the constant selection criterion that 
we planned to impose, since it had substantially 
higher total activity. We therefore neutrally 
evolved it for 16 generations in order to cre- 
ate P450s that were better equilibrated to the 
selection criterion. We evolved two parallel pop- 
ulations, which we named Rl and R2. The pro- 
cedure was exactly identical to that described 
below for the polymorphic evolution with the 
following exceptions: 

• Starting sequence: the starting sequence 
for the equilibration evolution was the 
21B3 sequence. 

• Population size: each of the two equili- 
bration evolution populations had a size 
of 174 sequences rather than the 435 used 
for the polymorphic evolution. 

• Selection criterion: the sequences were re- 
quired to have at least 75% of the total 
activity of the 21B3 P450. 



• Mutation rate: the mutation rate for the 
equilibration evolution was much higher 
than for the polymorphic evolution. The 
error-prone PGR protocol used 200 /iM 
manganese chloride (MnCl2), rather than 
the 25 fiM used for the polymorphic evo- 
lution. We estimate that this error-prone 
PGR protocol introduced ^ 4 nucleotide 
mutations per P450 gene at each genera- 
tion during the equilibration evolution. 

We performed 16 generations of equilibration 
evolution, and then randomly selected 23 func- 
tional mutants from each of the Rl and R2 
populations (Additional File [7]). We picked one 
of these mutants, Rl-11, for use as the parent 
for the neutral evolution experiments. 



Detailed protocol for evolution experiments 

We began with the Rl-11 P450 BM3 heme 
domain variant (sequence in Additional File 
[1]) cloned into the pCWori [25] plasmid with 
a 5' BamHl and 3' EcoRl site as described 
in [5]. The cloning primers were pCWoriJor (5'- 
GAAACAGGATCCATCGATGCTTAGGAGGTCAT- 
3' cmd pCWorLrev.clone (S'-GCTCATGTTTGACAGCTTATCATCG 
3'). We used error-prone PGR to generate 
mutants, taking great care to make the error- 
prone PGR protocol repeatable by using a rela- 
tively small number of thermal cycles. This was 
both to control the mutation rate by ensuring 
that the reaction did not saturate the reagents 
(which would cause the number of doublings to 
become sensitive to the initial template concen- 
tration), and to avoid the PGR-based recombi- 
nation events which can occur during with the 
last few thermal cycles of PGR reactions [36,37]. 
The PGR reactions were 100 ^1 in volume, and 
contained ~ 13 ng of plasmid template (cor- 
responding to Rr: 3 ng of template gene), 7 
mM magnesium chloride MgGl2, 1 x Applied 
BioSystems PGR Buffer II without MgGl2, 25 
pM. MnGl2, 0.5 pCWorLfor primer, 0.5 
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//M pCWori_rev primer, 200 /xM of dATP and 
dGTP, 500 Aim of dTTP and dCTP, and 5 units 
of Applied Biosystems AmpliTaq polymerase. 
The reactions were run on the BLOCK setting 
of a MJ Research PGR machine with a program 
of 95°C for 2 minutes, then 15 cycles of (95°C 
for 30 seconds, 57°C for 30 seconds, 72°C for 90 
seconds), and then cooling to 4°C. This proto- 
col yielded roughly 1-1.5 ^g of product gene (as 
quantified by gel electrophoresis versus a known 
standard), for a PGR efficiency of w 0.5. Se- 
quencing the unselected populations at the end 
of the experiment indicated that this protocol 
introduced an average of 1.4 it 0.2 nucleotide 
mutations, with the nucleotide error-spectrum 
shown in Table [TJ Because the number of PGR 
doublings is large compared the average muta- 
tion rate, the distribution of mutations among 
sequences should be well-described by the Pois- 
son distribution [38,39]. 

The mutant genes from the error-prone PGR 
were purified over a ZymoResearch DNA clean 
and concentrator column, and digested at 37°G 
with EcoRl and BamHl. The digested genes 
were then purified from an agarose gel with 
ZymoResearch DNA gel extraction columns, 
and ligated into pGWori plasmid that had been 
digested with BamHl and EcoRl and dephos- 
phorylated. The ligations were transformed 
into electro-competent catalase-free Escherichia 
coli [25] (the catalase is removed because it 
breaks down the hydrogen peroxide utilized by 
the P450 peroxygenase) , plated on Luria Broth 
(LB) plates containing 100 ^g/ml of ampicillin 
to select for the plasmid's antibiotic resistance 
marker, and grown at 37°G. Transformation of 
a control ligation reaction without any digested 
gene yielded at least 100-fold fewer colonies, 
indicating that the rate of plasmid self-ligation 
was less than one percent. 

Individual mutant colonies from the plates 



were picked into 96-well 2 ml deep-well plates 
containing 400 /il of LB supplemented with 100 
/ig/ml ampicillin. Each plate contained four 
parental control wells with cells carrying the 
parent Rl-11 gene, four null control wells with 
cells carrying the pGWori plasmid without a 
P450 gene, and a non-inoculated well to check 
for contamination. For the polymorphic pop- 
ulation, we picked five such plates with all 87 
other wells containing different mutants for a 
total population size of 5 x 87 = 435 mutants. 
For the 22 monomorphic populations (we be- 
gan with 24 populations but two had to be 
discarded due to contamination), we picked a 
single colony for growth and screening. For the 
unselected populations we picked a single colony 
for growth without screening. The LB deep-well 
plates were grown for 16-20 hours at 30°G, 210 
revolutions per minute (rpm), and 80% relative 
humidity in a Kuhner humidified shaker. To ex- 
press the P450 mutants, we prepared 2 ml deep 
well plates containing 400 ^1 per well of terrific 
broth (TB) supplemented with 200 iso- 
propyl /3-D-thiogalactoside (IPTG), 100 /^g/ml 
ampicillin, and 500 /iM of 5- aminolevulinic acid. 
We used a pipetting robot inoculated these TB 
plates with 100 //I from the LB plates. We 
stored the LB deep- well plates at 4°G, and 
grew the TB deep-well plates in the humidi- 
fied shaker at 30°G, 210 rpm, and 80% relative 
humidity for 22-24 hours. After this growth, 
the cells were harvested by centrifuging the 
TB deep- well plates at 4000 xg for 5 minutes 
and discarding the liquid. The cell pellets were 
flash- frozen in liquid nitrogen to aid in cell lysis. 



To lyse the cells for the assays, we resus- 
pended the ceh pellets in 300 ^ul of 100 mM [4- 
(2-hydroxyethyl)-l-piperazinepropanesulfonic 
acid] (EPPS) (pH 8.2) with 0.5 mg/ml lysozyme 
and 4 units /ml of deoxyribonuclease by pipet- 
ting 40 times with the pipetting robot. We then 
incubated the plates at 37° G for 30 minutes, 
again resuspended with the pipetting robot, and 
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put back at 37°C for an addition 30 minutes. 

We then pelleted the cell debris by centrifu- 
gation at 6000 xg for 5 minutes at 4°C. The 
pipetting robot was used to dispense 80 //I of 
the clarified lysate into 96-well microtiter plates 
(Rainin). We prepared a 6x stock of 1.5 mM 
12-pNCA in 36% dimethyl sulfoxide (DMSO) 
and the EPFS buffer (the 12-pNCA was stored 
in the DMSO solution and combined with the 
buffer immediately before use). We used a mul- 
tichannel pipette to add 20 fi\ of this substrate 
stock to each well of the microtiter plate. We 
briefly mixed the plates with "shake" setting 
of a 96-wcll plate spectrophotometer, and read 
an absorbance baseline at 398 nm. We then 
immediately added 20 /il of a freshly prepared 
solution of 24 mM hydrogen peroxide in the 
EPFS buffer to initiate the reaction, and mixed 
again. The final reaction conditions were there- 
fore the EPFS buffer with 6% DMSO, 4 mM 
hydrogen peroxide, and 250 fiM 12-pNCA. After 
40 minutes we quantified the amount of enzy- 
matic product by the increase in absorbance 
at 398 nm. This absorbance increase is due 
to the 4-nitrophenolate molecule released after 
the P450 hydroxylatcs the twelfth carbon of the 
12-pNCA molecule [23]. To score the mutants 
as functional or nonfunctional, we compared 
their gain in absorbance minus the median null 
control reading to that of the median parental 
control reading minus the median null control 
reading. All mutants that had at least 75% of 
the parental gain were scored as functional, all 
other mutants were scored as nonfunctional. 



We used the information from these assays 
to select the parents for the next generation. For 
the unselected population we did not require the 
mutants to be functional, so the selected mutant 
was used to start a 4 ml culture of LB with 100 
/xg/ml ampicillin, and the plasmid DNA was 
harvested with a mini-prep. This plasmid DNA 
was used as the template for the next round of 
error-prone PCR. Therefore, after the first gen- 



eration the four unselected replicates diverged 
into four separate error-prone PCR reactions. 
These unselected replicates were evolved for a 
total of twelve generations, and were sequenced 
at every third generation. 

For the polymorphic population, all mu- 
tants that were functional contributed an equal 
amount of plasmid DNA as template for the 
next generation. In order to do this, we col- 
lected 50 fA of the culture from the LB deep- well 
plate for each mutant that was scored as func- 
tional. All of these LB aliquots were pooled, 
and then the plasmid DNA was collected with 
a mini-prep. The pool of plasmid DNA was 
used as template for the next generation's error- 
prone PCR reactions. We performed 15 gener- 
ations of evolution for this polymorphic pop- 
ulation. Note that at each generation we are 
assaying 435 mutants as part of the evolution- 
ary procedure, so this provides information on 
mutational robustness. At every third gen- 
eration, we also selected a random sample of 
functional mutants for sequencing. After 15 
generations, we randomly selected 22 mutants 
for stability measurements and sequencing anal- 
ysis. The random selections were made from all 
functional mutants with the Python computer 
language random number generator. 

For the monomorphic populations, at each 
generation we assayed just a single mutant. If 
that mutant was nonfunctional, then at that 
generation the population stayed at its original 
sequence. In that case, for the next generation 
we simply picked a new mutant from the previ- 
ous generation's plate of transformed mutants. 
If the mutant we screened was functional, then 
that mutant represented the new population. 
We therefore grew a 4 ml LB culture with 100 
/ig/ml of ampicillin, and collected the plasmid 
DNA with a miniprep. That plasmid DNA was 
then used as the template for the next genera- 
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tion's error-prone PGR reaction. We thus had 
22 (actually 24 before 2 were contaminated) in- 
dependent monomorphic populations that were 
being evolved in parallel. Each was evolved for 
25 generations, and at the end of these 25 gen- 
erations we measured the stability of the final 
sequence of each population. Each time an as- 
sayed mutant was functional, we sequenced the 
new P450 gene. We also measured the average 
mutational robustness of the monomorphic pop- 
ulations at every fifth generation. To do this, 
we did a pooled mini-prep of equal volumes 
of LB cultures of all 22 replicates to obtain a 
equal mix of plasmid DNA. We then performed 
error-prone PGR on this mix, and assayed 435 
mutants to measure the fraction functional. Full 
neutral evolution data are in Additional FileO 



recombination. 



Test for recombination during error-prone 
PCR 

During the polymorphic population evolution, 
we performed error-prone PCR on a mix of 
different plasmids. It is common for PGR 
on mixed templates to lead to recombination 
events during the reaction [36, 37]. We at- 
tempted to reduce this recombination by using 
a small number of thermal cycles. However, in 
order to test for recombination, we analyzed 
the sequences of the final 22 selected members 
of the polymorphic population. There are a 
variety of statistical tests to detect recombi- 
nation in a set of sequences. A comparison 
of these tests by Posada [40] found that the 
Max-Ghi^ method developed by John Maynard 
Smith [41] performs well. A publicly avail- 
able implementation of this method [42] is at 
http : / / www . lifesci . Sussex . ac . uk/ C SE/test / maxchi , 
We used this implementation to analyze the 22 
final polymorphic sequences, and the resulting 
P-value was 0.29 after 100 random permuta- 
tions, indicating that there is not significant 



Measurement of P450 stabilities 

We measured the stabilities to both irreversible 
thermal and irreversible urea denaturation 
of the final (generation 25) member of each 
monomorphic population, as well as of the 22 
randomly selected members of the polymorphic 
population. As discussed in the Supplemen- 
tary Information of [5], cytochrome P450 BM3 
heme domains (and indeed most P450s) dena- 
ture irreversibly, forcing us to use resistance 
to irreversible denaturation to quantify protein 
stability. The first stability measure is the Tsq, 
defined as the temperature at which half of the 
protein is denatured after a 10 minute incu- 
bation. The second stability measure is the 
[urea] 50, defined as the urea concentration at 
which half of the protein denatures after a 4 
hour room-temperature incubation. Each set 
of measurements (those of T50 and [urea] 50) 
was performed on all of the mutants in the 
same day, and each mutant was treated identi- 
cally. Therefore, it is possible to make accurate 
comparisons of the relative values of the mea- 
surements within the data set. However, the 
absolute values of the T50 and [urea] 50 values 
may be less accurate. Therefore, care should be 
taken in comparing the absolute value of these 
measurements to those of other studies (such 
as [5]). 



Both the T50 and [urea] 50 measurements 
were performed in clarified cell lysate. The 
protein was expressed using catalase-free E. 
coll [25] containing the encoding gene on the 
]JE[3!G inducible pGWori [25] plasmid. We used 
freshly streaked cells to inoculate 2 ml cultures 
of LB supplemented with 100 /Ug/ml of ampi- 
cillin, and grew these starter cultures overnight 
with shaking at 37° G. We then used 0.5 ml 
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from these starter cultures to inoculate 1 L 
flasks containing 200 ml of TB supplemented 
with 100 /ig/ml of ampicillin. The TB cultures 
were grown at 30°C and 210 rpm until they 
reached an optical density at 600 nm of ?a0.9, 
at which point IPTG and (5-aminolevulinic acid 
were added to a final concentration of 0.5 mM 
each. The cultures were grown for an additional 
19 hours, then the cells were harvested by pel- 
letting 50 ml aliquots at 5,500 g and 4°C for 10 
min, and stored at -20°C. To obtain clarified 
lysate, each pellet was resuspended in 8 ml of 
100 mM EPFS (pH 8.2) and lysed by sonica- 
tion, while being kept on ice. The cell debris 
was pelleted by centrifugation at 8,000 g and 
4°C for 10 minutes, and the clarified lysate was 
decanted and kept on ice. 



For the T50 measurements, 125 fil of clar- 
ified lysate from a single mutant was added 
to all 12 wells in a row of a 96-well hard-shell 
thin- wall microplate (MJ Research). The plate 
was heated for 10 minutes using the gradient 
method of an Eppendorf Mastercycler gradient 
FCR machine, with the gradient set at either 
33°C-45°C or 46°C-58°C (each mutant was ex- 
posed to both of these gradients), the machine 
on the BLOCK setting, and the heated lid set to 
75°C with the lid WAIT option. The plate was 
then cooled to 4°C, removed from the FCR ma- 
chine, and centrifuged at 5,500 g and 4°C for 5 
minutes to pellet any debris. A pipetting robot 
was used to dispense 80 fil of the supernatent 
into a 96-well microtiter plate (Rainin) , and the 
amount of remaining properly folded F450 was 
quantified from the carbon monoxide difference 
spectrum as described below. The T^q values 
were determined by fitting sigmoidal curves the 
percent of remaining protein as shown in Addi- 
tional File [3l Our ability to accurately compare 
T50 values within the data set requires that each 
well in a given column of the gradient FCR ma- 
chine be at the same temperature. We used 
a thermocouple to measure the temperature of 



the wells with the machine lid open, and con- 
firmed that the wells were within a few tenths 
of a degree of the same temperature. Further 
evidence for the consistency of our T50 val- 
ues comes from the fact that two independent 
measurements of the T50 for our Rl-11 par- 
ent yielded values that differed by only 0.1°C. 
However, the absolute values of the measured 
temperatures are less accurate. Thermocouple 
measurements indicated that, with the machine 
lid open, the wells were ^ 1°C cooler than the 
indicated temperature. We were unable to as- 
certain the temperatures with the heated lid 
closed, but based on comparisons water bath 
measurements, the temperatures with the lid 
closed slightly exceeded the indicated tempera- 
tures. 



For the [ureajso measurements, 125 /ul of the 
clarified lysate from a single mutant was added 
to all 12 wells in a row of a 96-well microtiter 
plate. A pipetting robot was then used to add 
and mix 125 fil of a 2X solution of urea in 100 
mM EFFS (pH 8.2) so that each subsequent 
column had a higher concentration of urea, and 
so that the final urea concentrations were those 
shown in Additional File [H The plates were 
left on the bench at room temperature for 4 
hours, and the amount of remaining properly 
folded F450 was quantified from the carbon 
monoxide difference spectrum as described be- 
low. The [urea] 50 values were determined by 
fitting sigmoidal curves to the percent of re- 
maining protein. Evidence for the consistency 
of the [urea] 50 measurements comes from the 
fact that two independent measurements of the 
[urea] 50 for our Rl-11 parent yielded values that 
differed by only 0.01 M. In addition, the [urea]5o 
and T50 values are highly correlated (Additional 
File [5|) , indicating that they provide consistent 
measures of stability. 

For both the r5o and [urea]5o measurements. 
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the folded P450 was quantified from the car- 
bon monoxide difference spectrum [43]. The 
microtiter plates containing the P450 samples 
were first used to read blank spectra at 450 and 
490 nm using a Tecan Satire 2 plate reader. The 
plates were then incubated for 10 minutes in an 
airtight oven with carbon monoxide. The plates 
were removed form the oven and 10 /xl of 0.1 M 



sodium hydrosulfite in 1.3 M potassium phos- 
phate (pH 8.0) was immediately added to each 
well. After 5-10 minutes, spectra were again 
read at 450 and 490 nm. The amount of P450 
is proportional to the increase in the signal at 
450 nm after this procedure minus the change 
in the signal at 490 nm. 
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A.l Mathematical background 

The first purpose of this appendix is to provide mathematical equations that describe the experi- 
ments. The second is to show how four measurements from the experiments can be used to calculate 
two quantities that describe the topology of the underlying protein neutral network. We will derive 
two equations for both quantitites, each in terms of a different measurement. The fact that the four 
equations will be seen to yield consistent results provides evidence for the accuracy of the following 
calculations. 



Our calculations are based on a view of neutral protein evolution as a process constrained by a 
stability threshold, a view that we originally introduced to explain experimental protein mutagen- 
esis results [3]. The calculations closely parallel our earlier work [22], which is in turn based on a 
general theoretical treatment of evolution on neutral networks by van Nimwegen and coworkers [11]. 
These calculations will probably be most thoroughly understood by first reading those works. The 
primary difference between the current calculations and [22] is that previously we assumed that 
the per generation per protein mutation rate fi was <C 1, so that at each generation a protein was 
either unmutated (with probability 1 — fJ.) or experienced a single mutation (with probability /i). In 
contrast, here we allow the mutation rate to be arbitrarily large, so that a protein may experience 
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multiple mutations in a single generation (in this sense the calculations resemble the generalization 
by Wilke [18] of [11]). Specifically, let fm be the probability that a protein experiences m mutations 
in a single generation. Here we derive results for arbitrary and then approximations relevant to 
the form of fm in the experiments. In the limiting case of small mutation rate (where /o = 1 — 
/i = n, and fm = fov m > 1), the calculations here reduce to those in [22]. Proteins evolving 
in nature typically experience very low mutation rates, so [22] probably offers the best description 
of natural protein evolution. The calculations presented here are designed to specifically treat the 
evolutionary dynamics of the experiments. 

A protein's thermodynamic stability is described by its free energy of folding, AG/, with more 
negative values indicating more stable proteins. As described in several previous papers [3,5,22], 
we assume that selection requires a protein to fold with some minimal stability AG™™, so that a 
protein adequately folds if and only if AG/ < AG™™. The amount of extra stability a protein 
possesses relative to the stability threshold is given by AG/'*'"^ = AG/ — AG™™; note that all 
folded proteins will have AG^^^^ < 0. Wc further assume that as long as AG/^*'^ < 0, selection is 
indifferent to the exact amount of extra stability that a protein possesses (see [22] for a discussion 
of the limitations of this assumption). We conceptually divide the continuous variable of protein 
stability into small discrete bins of width b. Specifically, a protein is in bin i if it has AGJ^^^ be- 
tween (1 — i)b and —ib, where i = 1, 2, . . .. Mutating a protein changes its stability by an amount 
A AG (defined as the stability of the mutant protein minus the stability of the initial protein) , and 
so may move it to a new stability bin. In [22], we defined a matrix W with elements Wij giving 
the transition probabilities that a single mutation changes a protein's stability from bin j to bin i. 
We noted that W could be computed from the distribution of AAG values for all single mutations, 
and argued that W remains fairly constant during neutral evolution since the distribution of AAG 
values remains relatively unchanged. However, we emphasize that (as discussed in detail in [22]) 
the constancy of the AAG distribution remains an assumption, albeit one that has now been shown 
to be quite accurate for lattice proteins [3,22,44] and provide a consistent theoretical explanation 
for a growing body of experimental results (the current work as well as [3]). 

Since we are allowing for larger mutation rates, and we must consider the possibility that a 
protein's stability might change due to multiple mutations at a single generation. Therefore, we 
make a more general definition of Wij^m as the probability that m random mutations to a protein 
in stability bin j change its stability to bin i, and let Wm be the matrix with elements Wij,m. Note 
that Wm only describes mutations that cause transitions from one folded protein to another, since 
the stability bins z = 1, 2, ... all correspond to folded proteins. As before [22], we assume that Wm 
is roughly constant during evolution, meaning that the distribution of AAG values for multiple 
mutations is roughly constant during neutral evolution. Note that if m = 1, then Wm is just the 
matrix W that can be computed from the distribution of single-mutant AAG values [22]. We will 
now use the matrices Wm to calculate the following characteristics of a population that has evolved 
to equilibrium: the distribution of stabilities, the average number of mutations {m)T accumulated 
after T generations, and the average fraction {T) of stably folded proteins in the population. We 
then introduce a few approximations (that should be quite accurate for the experimental work in 
this paper) that greatly simplify these calculations. Finally, we relate the calculations to properties 
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of the underlying protein neutral network. 



As described generally by van Nimwegen and coworkers [11], the evolutionary dynamics depend 
on whether the evolving population tends to be monomorphic or highly polymorphic. When the per 
sequence per generation mutation rate /U is <C 1, whether the population is mostly monomorphic or 
highly polymorphic is determined by the product of the population size A'^ and /i: when N <^ 1 
the population is mostly monomorphic, and when Nfi ^ 1 the population is highly polymorphic. 
However, with multiple mutations per generation, Nfi is no longer an appropriate parameter to 
distinguish between mono- and polymorphism, since if the population size is sufficiently small the 
population can still be monomorphic even if there are multiple mutations per generation. Specifi- 
cally, in one set of experiments we constrained the population to be monomorphic (by maintaining 
a population size of one), but still allowed the single protein in this population to experience more 
than one mutation at a generation. So we instead denote the populations as either monomorphic or 
polymorphic. We indicate quantities calculated for the monomorphic population by the subscript 
M (i.e. {T)m) and those calculated for the polymorphic population by the subscript P (i.e. {!F)p). 



A. 2 Monomorphic limit 

In the limit of a completely monomorphic population, all of the proteins are in a single stability 
bin. Let pi (t) be the probability that the population is in stability bin i at time t, and let p (t) 
be the column vector with elements pi (t). At each generation there is a probability /o that there 



oo 



is no mutation that becomes fixed in the population, a probability of Yl fmYl ^ij,mPj that the 

m=l j 

population experiences a mutational event (which could be a single mutation or several simultaneous 



oo 



mutations) that moves it into bin i, and a probability Yl fmPi Yl ^ji,m that the population is in 

m=l j 

bin i and experiences one or more mutations that move it to another bin of stably folded proteins. 
Define i/j^m = Y2 ^ji,m to be the fraction of m-mutants of a protein in bin i that still fold, and let 

j 

Vm be the matrix with diagonal elements given by Va^rn = '^i,m and all other elements zero. The 
time evolution of p is 



P (t + 1) 



Pit) (1) 



I + 5^ /m (W„, - V„,) 

m=l 

where I is the identity matrix. Note that mutations that destabilize a protein beyond the stabil- 
ity threshold are immediately lost to natural selection, and so leave the population in its original 
stability bin. This describes the experiments for the monomorphic populations, where we retain 
the parental sequence if the single mutant we generate is nonfunctional. Equation [T] corresponds to 
Equation (1) of [22], and the blind ant random walk described by van Nimwegen and coworkers [11]. 



Equation [T] describes a Markov process with a non-negative, irreducible, and acyclic transition 
matrix, and so p approaches a unique stationary distribution (equilibrium value) of pm given by 
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the eigenvector equation 



Pm 



m=l 



PM- 



(2) 



Once p has reached equihbrium, the average fraction of proteins that still stably fold at each 
generation is 



m I Pm 

where e = (l,...,l)is the unit row vector. 



(3) 



m=l 



To calculate {m)T,M, the average number of mutations accumulated after T generations once the 
population has equilibrated, we note that at each generation there is a probability of fmPj ^ij,m 

i 

that a randomly chosen protein is in bin j, experiences m mutations, and still stably folds. The 
average number of mutations accumulated in a single generation is simply the average of m weighted 
over this probability. So summing over all values of m and j, we see that 

oo 

{m)T,M = Te ^ m/^WmPM- (4) 

m=0 

This equation corresponds to Equation (6) of [22] , which was derived using an embedded Markov 
process formalism. Here wc have foregone this formalism for the more intuitive argument presented 
above, since we do not attempt to calculate higher moments of the number of mutations. 



A. 3 Polymorphic limit 

In the limit when the population is highly polymorphic, at each generation there are sequences in 
many different stability bins. In this case, we describe the distribution of stabilities by the column 
vector x(t), with element Xi(t) giving the fraction of proteins in stability bin i at time t. At 
generation t, the fraction of mutants that continue to fold is 

(^)t = e ( /ol + ) X (t) . (5) 

V m=l / 

Therefore, in order to maintain a constant population size, each remaining protein must produce 
an average of at = {^)t~^ offspring. The population therefore evolves according to 

x(t + l) = aJ/oI+f;/^W™)x(t). (6) 

V m=l / 

After the population evolves for a sufficiently long period of time, x will approach an equilibrium 
value of xp. At this equilibrium, the average fraction of mutants that fold at each generation is 

{J^)p = elfoI+f2 /m ) xp , (7) 

\ m=l J 
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and the equilibrium reproduction rate is q = {J^)p ^ ■ Therefore, 



xp = a /ol + ^ /^W^ xp. (8) 



Equations[7|and[8]can be combined to show that xp and {J-) p can be calculated from the eigenvector 
equation 

oo 

{{T)p - fo) xp = /mW^xp, (9) 

m=l 

oo 

with {{J-)p — fo) the principal eigenvalue of the nonnegative and irreducible matrix fm^^m- 

m=l 

Equation [9] corresponds to Equation (14) of [22], Equation (6) of the work by van Nimwegen and 
coworkers [11], and Equation (13) of the work by Wilke [18]. 



We now calculate {m)T,p, the average number of mutations accumulated in T generations after 
the population has equilibrated. At equilibrium, there is a probability of fmXj Yl ^ij,m that a 

i 

protein is in bin j, experiences m mutations, and still stably folds. Subsequently, all of these folded 
proteins produce an average of a offspring. The average number of mutations accumulated in a 
single generation is simply the average of m weighted over this probability, and then multiplied by 
the average reproduction rate. So summing over all values of m and j, we obtain 

oo J, oo 

{m)T,p = aTe ^ m/^W^xp = "t^:,— e ^ m/^WmXp. (10) 

m=0 ^ m=0 

This equation is the counterpart of Equation (18) of [22], where we have again foregone the em- 
bedded Markov process formalism for a more intuitive derivation. 



A. 4 Approximations for polymorphic limit 

We can dramatically simplify the results from the previous sections with several reasonable ap- 
proximations. The first approximation is that the AAG values for random mutations are roughly 
additive, and is supported by a number of experimental studies of the thermodynamic effects of 
mutations [45-47]. We have previously shown that this approximation can be used to accurately 
describe experimental protein mutagenesis data with a simple stability threshold model [3] . Under 
this approximation, the distribution of net AAG values for multiple mutations can be computed 
from the distribution of AAG values for single mutations by performing convolutions of the single- 
mutation A AG distribution [3] , meaning that for arbitrary m can be computed solely from the 
distribution of AAG values for single mutations. However, to simplify the equations from previous 
sections, we need to express Wm for arbitrary m only in terms of W (recall that W = Wi). Since 
W only contains information about stability transitions from folded proteins to other folded pro- 
teins, if we make the second approximation that a protein that is destabilized beyond the minimal 



16 



stability threshold by one mutation is not re-stabilized to a folded protein by a subsequent muta- 
tion, then Wjn = W". This approximation that unfolded proteins are not re-stabilized should be 
quite accurate since stabilizing mutations tend to be relatively rare and small in magnitude [48-51] 
(this is the underlying idea behind the Markov chain approximation that was shown to be highly 
accurate for lattice proteins [44]). To summarize, if A AG values are roughly additive and stabilizing 
mutations are rare, we have the approximation 

Wn, « W™. (11) 



Simplifying the equations of the previous sections also requires assigning a specific functional 
form to fm, the probability that a sequence undergoes m mutations. Here we assume that mutations 
are Poisson distributed among sequences, so that 



fm = (12) 

m! 

oo 

where /U = ^ i^fm is the average number of mutations per protein per generation. When the mu- 

m=0 

tations are introduced by error-prone PGR, the Poisson distribution is an excellent approximation 
to the true theoretical distribution of mutations created by error-prone PGR [38, 39] provided that 
^ is much less than the number of PGR doublings, as is the case in all of the experiments in the 
current work. 



We now use the approximations of Equations [TT] and [T2] to simplify the results given above for 
the highly polymorphic limit. We begin by using these approximations to rewrite Equation [9] as 



oo 



((.F)p - e-^) xp = e-'^^^W-xp. (13) 



m! 

m=l 



This equation makes clear that xp is the principal eigenvector of the matrix ^ ^W™", therefore 

m=l 

Xp must also be the principal eigenvector of W. Now in our earlier work [22], we defined the 
principal eigenvector of W as Xqo, called the corresponding eigenvalue (z^)oo) and showed that this 
eigenvalue is shown the average fraction of single mutations that are neutral in a population that 
is evolving with Njj, ^ 1 and ^ 1. Therefore, with the approximation of Equation [TTl xp and 
Xoo are equal, and are both defined by the same eigenvector equation, 

(zy)ooXp = Wxp = WXoo = (i^)ooXoo- (14) 

Combining Equations [13] and [14] we have, 

(y-)pxp — e 2^ — — — 



m,=0 



oo 



xp (15) 
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Equation [15] can be solved to yield 



(.)oo = l + i^. (16) 



Similarly, we can simplify Equation IIOI 



rp oo 



^-^ ml 

m=l 

^ ml 

m=l 

^ — ' m 



ml 

m=0 

= Tii{u)^. (17) 
Solving this equation for {v)oo yields 

M« = (18) 



A. 5 Approximations for monomorphic limit 

We now simplify the equations for the monomorphic limit. This requires several further approxima- 
tions. We begin by approximating that the stability probability distribution pM given by Equation 
[2] by the distribution po defined in [22] as satisfying 

= (W-V)po. (19) 

The basic rationale behind approximating pM with po is that Equation [2] can be viewed as a per- 
turbation to Equation [19] [52]. Essentially, po is an eigenvector of the matrix W — V while pM is 

OO rn — 1 

the corresponding eigenvector of the matrix W — V + ^ - — j— (W" — Vm). The latter matrix 

m=2 ™" 

OO rn — l 

can be viewed as a perturbation to the first, since the sum ^ — j— (W" — Vm) is small. This 

m=2 

smallness is due to the fact that W™' tends to zero with large m, causing Vm to tend towards the 
identity matrix. In addition, the f/^ /m\ terms tend to zero with large m. Therefore, the terms in 
the summation are all simply either a perturbation to W — V or involve subtracting terms that 
are fractions of the identity matrix. The perturbations lead to bounded changes in the eigenvec- 
tors [52], while the identity matrix terms do not change the eigenvectors. Below we give a more 
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rigorous justification of the assumption that pM is approximately equal to po- 



We need one additional approximation to make further progress. Both Equations [3] and U] con- 
tain terms of the form WmPoj and even if we use Equation [11] to rewrite these terms as W™po, 
there are no further clear simplifications. However, any probability vector that is multiplied re- 
peatedly by W and normalized will eventually converge to x^o = xp (since this is the principal 
eigenvector of W). We make the approximation that this convergence is sufficiently rapid to be 
essentially complete after a single multiplication. This approximation is supported by both protein 
mutagenesis studies [3,26,27] that indicate that proteins rapidly converge to an exponential decline 
in the fraction folded (indicating the stability distribution has equilibrated, as discussed below, and 
by lattice protein studies showing the same [3,44]. Therefore, we make the approximation that 
eW'^po = (i^)oeW'"^"'^Xoo = {v)o{v)ac!^~^ where {v)o = eWpo has the same definition as in [22], 
where it was defined as the average fraction of functional single mutants of a population evolving 
with ^ < 1 and N^x < 1. 

We use these approximations to simplify Equation [3] as 





oo 



'o 



m=l 




oo 



m—l 




(20) 



Solving this equation for (j^)o, we find 




(21) 
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We now use the approximations to simplify Equation [5] as 

CO 

{m)T,M = Te'^ mfmW^pM 



Solving this equation for {i')o yields 



m=0 

oo 



^-^ ml 

m=l 

°° m 

Z / JT-l 1 



ml 

171=1 



/iT(i/)oe^(^''^~~^). (22) 



To recap, we now have equations to calculate (i^)oo and {u)o from experimentally measurable 
quantities. Equations [T6l and [TSl allow us to calculate {v)oo from (J^)p and {m)T,p, respectively. 
Given this calculated value of {v)oo, Equations [21] and [23] then allow us to calculate {u)^ from {J-)m 
and {m)T,M, respectively. The fact that we have two equations each for (i^)oo and (z/)^ allows us 
to assess the self-consistency of the approach. 



A. 6 Interpretation in terms of neutral networks 

Throughout the preceding calculations, we have referred to (z^)oo and (i^)o as we defined them 
in [22]: namely, as the average neutrality of protein populations evolving with /i <C 1 and Nfj, 
either 3> 1 or <C 1, respectively. However, van Nimwegen and coworkers [11] have shown that 
they can also be interpreted in terms of the underlying neutral network. In the experiments we 
make mutations at the nucleotide (rather than amino acid) level, so each point in our sequence 
space corresponds to a different gene. Every gene that yields an amount of protein sufficient to 
hydroxylate the twelfth carbon of 12-p-nitrophenoxydodecanoic acid with at least 75% of the total 
activity conferred by the original Rl-11 parent gene represents a node on this neutral network. 
We note that in the experiments (and also usually in natural evolution), the edges on the neutral 
network are not all completely equivalent or fully undirected, since some mutations are more likely 
to occur than others (for example, error-prone PGR with Taq polymerase is more likely to cause an 
A^G mutation than an A— s-G mutation). In the analysis that follows, we ignore this complication 
and assume all neutral network edges are equivalent. 
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In an extremely insightful analysis, van Nimwegen and coworkers [11] have shown that impor- 
tant characteristics of a neutral network can be inferred from evolutionary quantities. Specifically, 
they have shown that if a population is evolving with <C 1 and N (j. ^ 1, then the average 
neutrality (which we have denoted by (z^)oo) is equal to the principal eigenvalue of the adjacency 
matrix of the neutral network, normalized by the network coordination number (number of possi- 
ble connections per node). In addition, they pointed out that a population evolving with /x <C 1 
and N ^ <C 1 moves like a blind ant random walk, meaning that the average neutrality (which we 
have denoted by {v)o) is equal to the average connectivity of a neutral network node divided by 
the network coordination number. In our P450 experiments, we have measured the values needed 
to estimate (i^)oo and {u)o using Equations [161 HI EH and [23l Using the final values listed in 
Table [21 {!F)p = 0.50 and {!F)m = 0.39. Taking the final nucleotide mutation values from Table 
m {m)T,p/T = 0.69 and {m)T,M/T = 0.31. The average mutation rate, computed from the un- 
selected population, is /i = 1.40. So using Equation 1161 ( 

i^) — 0.53, while using Equation 1181 
{t^)oo = 0.49. The consistency of these two values supports the idea that the calculations above 
accurately describe the evolutionary process. Taking the average value of these two measurement 
as {1^)00 = 0.51, we can then use Equations [2T] and [23] to calculate (i/)o. We calculate values of 
0.28 and 0.43, respectively. These estimates differ by more than those for (z^)oo) perhaps because 
additional approximations have gone into the derivation of the relevant equations (in addition, we 
have made no attempt to carry out the rather complex propagation of the sampling errors of Table 
[2|). However, the values are still in a similar range. Taking the average of these two values, we 
estimate that {i')o = 0.35. So overall, we predict that each functional P450 gene should have an 
average fraction of 0.35 of its sequence nearest neighbors also encoding a functional gene, for an 
average of about 1,500 neighbor genes. We predict that the principal eigenvalue of the neutral 
network adjacency matrix is 0.51 x3L. The fact that principal eigenvalue exceeds the average 
connectivity indicates that the neutral network is not a regular graph, but instead has some nodes 
more highly connected than others. 

The value for (z^)oo calculated above can also be related to measurements from protein mu- 
tagenesis experiments. A number of studies [3, 26, 27] have observed that the probability that a 
protein remains functional after m mutations falls off exponentially with the number of mutations. 
In fact, the decline is not always exponential for the first few mutations if the starting protein has 
especially high or low stability [3] or activity [53], but will still converge to this exponential form 
after a few mutations [3,44,54]. The stability threshold model can be used to relate this decline to 
(z^)oo, as is done indirectly in the Markov chain approximation of [44]. Here we make that connec- 
tion explicit. The initial protein has a stability that falls into some stability bin i. Therefore, its 
stability can be described by the column vector yo, which has element i equal to one and all other 
elements equal to zero. Now imagine constructing all single mutants of this protein. The fraction 
of these single mutants that still fold is just eWyo, and the distribution of stabilities among the 
single mutants is yi = Wyo (note that the elements of yi no longer sum to one). Similarly, after 
m mutations, the fraction of mutants that still fold is eWmyo, and the distribution of stabilities 
among the m-mutants is ym = Wmyo. With the approximation of Equation \TT\ = W^yo. 
This makes it clear that ym will converge to a vector proportional to Xqo, the principal eigenvector 
of W. Once this convergence is complete, each new mutation simply reduces the fraction of mu- 
tants that fold by a factor of {1^)00, the principal eigenvalue of W (and the spectral radius of the 
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neutral network normalized by the coordination number). Therefore, what we have called (z^)oo in 
the present work and [22] is equal to what is called x in [27], q in [26], and (u) in [3]. The major 
difficulty that is faced in extracting (z^)oo by the method of those three studies [3,26,27] is that it 
is not possible to directly assay mutants with m mutations, but instead only possible to assay a set 
of mutants with a distribution of m. All three studies use different (and valid) methods to account 
for this distribution, but this accounting is still difficult because most of the functional mutants 
come from the low m end of the distribution. This makes it hard to get accurate values for the 
fraction functional after large numbers of mutations, since most of the functional mutants in the 
set come from sequences with few mutations. For this reason, we believe the current method of 
measuring (j^)oo is more accurate. A second caution about comparing values of (i^)oo from different 
studies is that its value depends on the nucleotide error-spectrum of the experiment, since differ- 
ent mutagenesis methods create different distributions of nucleotide and amino acid mutation types. 

We also briefly mention how we arrived at an estimate of (i^)oo for 3-methyladenine DNA 
glycosylase from the data of [27]. This paper reports that a fraction x = 0.34 of amino acid 
mutations inactivate the protein. We would like to determine the fraction (z^)oo of nucleotide mu- 
tations that do not inactivate the protein. Roughly 75% of random mutations to a gene will be 
synonymous. Therefore, m amino acid mutations should cause about 4m/3 nucleotide mutations. 
The study of [27] measures that after m mutations, a fraction (1 — x)^ of the mutants are func- 
tional. That means that (1^)00^'"''^ fraction should be functional. Equating these expressions yields 
{1^)00 = exp (I log (1 — x)). So using x = 0.34, we arrive at {i')oo = 0.73. 



A. 7 Detailed justification for approximating pM by po 

Here we provide a detailed justification for the approximation that pM is about equal to po. In the 
monomorphic limit, the time evolution of p is given by Equation [H and the stationary distribution 
Pm is given by Equation [2j We assume the approximations of Equations [11] and [12] and show that 
we can approximate pm by po, where po is given by Equation [19] To justify this approximation, 
we insert po into the right hand side of Equation [T] and ask to what extent po is left unaltered by 
the dynamics. If po is found to be stationary to good approximation then, by uniqueness of the 
stationary distribution of an ergodic process, po would be a good approximation to pM- 

We therefore suppose that at some time t the distribution is given by po and compute, using 
Equation [T] the change in po after one generation. The new distribution at time t + 1 is given by 




(24) 



Using (V — W) Po = 0, and taking components of the above equation, we obtain 



00 




(25) 



m=2 
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Thus Po would be an approximately stationary distribution of the dynamics if | X] fm [(W"^ — Vm) Po]J <S 

m=2 

PQi. We now proceed to show that this will be the case in most situations of interest by deriving 
upper and lower bounds on the second term of the right hand side of Equation [25l 



Consider first the term (W'"po)j, which can be written as 



POk 

T7 



where we have used Wpo = Vpo in the second equality. We now note that Uk < fmax for all 
k, where i^max is the maximum neutrality, maximized over all bins. This leads to the successive 
inequalities 

fcl , . . .,km — l 
kl,...,km-2 

< ^max Y ^iki^kiki ' " " M^fc™_3fc„_2P0fc™_2 

fcl,...,fcm-2 

<<ax'E^*'^i^Ofci, (27) 



yielding the upper bound 



W™Po). < Ca"x'^^.PO^. (28) 



/I — max 

In an identical manner, we obtain the lower bound 

(W™po), > u^X\po^, (29) 

where Umin is the smallest neutrality, minimized over all bins. Note that both inequalities above 
become exact equalities when all bins have the same neutrality v, which could be interpreted as 
either u^m or z^max- 



Having obtained inequality constraints on (W"*po)j, we now consider the term (VinPo)j! which 
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can be written as 

(VmPo)j = Pom,m 



= Poi Yl ^i^i ^^1^2 ■ ■ ■ ^fc— 1« 

= POi Yl ^fci^fcifca 

^1 vi^rn— 1 

hi , . . .,km--i 



which yields an identical upper bound to that on (W™po)j, namely, 



li - ""max '^iFUi-, 

and similarly 



(V-po),, < i^^^-^Poi, (31) 



(V-po), > ^.Poi. (32) 

It should again be noted that both the above inequalities become exact equalities when all bins 
have a common neutrality z^. 



We are now in a position to estimate bounds on the magnitude of the second term of Equation 
[25l Using the four inequalities of Equations [28l [29l EU and [32] above, we have 

- K^.' - C;') ^^P0^ < [(W- - V^) po], < - C;') ^^P0^, (33) 

or equivalently, 

|[(W--V^)poy < «ax' - C;') ^^Po., (34) 

where the inequality above becomes an exact equality when all bins have the same neutrality. How- 
ever, in this limit, the right hand side of the above equation vanishes, and therefore the second 
term of Equation [25] is identically zero in this case, giving the result that pM is exactly equal to 
Po when all bins have the same neutrality, even if fj, is arbitrarily large. 



We now carry out the sum over m to obtain an upper bound on the second term of Equation [25l 
in the more general and realistic case of unequal neutrality bins. Using Equation [M] and the specific 
Poisson form of fm, we obtain an upper bound on the fractional change in poi in one generation: 



Piit + 1) 



POi 



POi 



,m—l 



1,11 

^ — ^ m! 

m=2 



,m—l 



) 



(35) 
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The above bound vanishes for smah /i, is an increasing function of fmax — z^min; and is typically 
much smaller than 1. An extreme estimate of the size of the fractional change can be made when 
t'max = 1 and fmin = 0. In this case, using fi = 1.4 (the value in our experiments), the above 
inequality simplifies to 

< i/i (l - e-f" - /xe"^) ~ OAli^i. (36) 

Noting that i^j < 1, the fractional change in poi is therefore reasonably controlled even in the most 
extreme case. For realistic situations, the fractional change in pQi is expected to be much lower, 
thus justifying the use of po as the stationary distribution of the dynamics of Equation [TJ 



Pi{t + 1) -poi 
POi 
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Figure 1 - Theoretical views of the evolution of protein mutational robustness and stability. 

(A) Theory predicts that a mostly monomorphic population is equally likely to occupy any node of its neutral 
network, while a highly polymorphic population will prefer more connected nodes [11]. Node sizes are drawn 
proportional to the occupation probabilities. (B) Proteins evolving in a highly polymorphic population 
are predicted to be more stable than their counterparts in a mostly monomorphic population [22]. The 
histograms illustrate the distributions of stabilities for the two cases. The increased stability is a biophysical 
manifestation of excess mutational robustness, since more stable proteins are more mutationally robust [3-5]. 
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Figure 2 - Outline of the neutral evolution experimental procedure. 

For the polymorphic population, error-prone PGR was used to generate mutant P450 genes. These genes 
were hgated into a plasmid and transformed into E. coli. Individual mutants (435) were picked, expressed 
in E. coli, and assayed for enzymatic activity. All mutants that met the selection criterion contributed an 
equal amount of plasmid DNA as template for the next generation of error-prone PGR. The monomorphic 
populations were treated similarly, except only a single mutant was assayed at each generation. If this mutant 
met the selection criterion then it became the template for the next generation of error-prone PGR; otherwise 
at the next generation another colony was picked from the same plate. In the unselected populations a 
single mutant was picked and used as the template for the next generation of error-prone PGR. 
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Figure 3 - Accumulation of nucleotide {{rrint)) and nonsynonymous {{niaa)) mutations in the 
experimentally evolved P450 populations. 

For the unselected and monomorphic populations, numbers are the average over all replicates at the indicated 
generation; for the polymorphic population they arc from a random sample, with sampling standard error 
shown. 
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Figure 4 - The polymorphic population neutrally evolved a higher average mutational robustness 
than the monomorphic populations. 

The fraction functional was determined by assaying 435 mutants (average of 1.5 nucleotide mutations per 
gene). Error bars show binomial standard error. For the monomorphic population, numbers are the average 
over all replicates. 
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Figure 5 - The more mutationally robust proteins are more stable. 

The P450s from the polymorphic population neutrally evolved higher stability and expression levels than 
their counterparts from the monomorphic populations. The histograms show the distributions for the final 
protein from all monomorphic replicates and for the same number of randomly chosen proteins from the final 
polymorphic population. The plots show (left to right) the temperature at which half the protein irreversibly 
denatured after 10 minutes (Iso), the urea concentration at which half the protein denatured after 4 hours 
([urea] 50), and the expression level relative to that of the original parental P450. The means are significantly 
different, with unequal variance t-test P-values of 0.02, 0.005, and 0.04, respectively. 
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Tables 

Table 1 - Error-prone PCR nucleotide mutation spectrum. 

Spectrum of nucleotide mutations introduced by the error-prone PCR procedure used in the neutral evolution 
experiments. The spectrum was determined by sequencing the four final (generation 12) sequences from the 
unselected population, since in these sequences the mutations accumulate without constraint. As has been 
previously noted for error-prone PCR with Taq polymerase [3,5,26], the nucleotide error spectrum is biased 
towards certain types of mutations. 



Total nucleotide mutations 


67 


% synonymous mutations 


25 


Mutation types (%) 




A ^T, T ^A 


19.4 


A ^C, T 


1.5 


A ^G, T ^C 


64.2 


G ^A, C ^T 


4.5 


G ^C, C ^G 


0.0 


G ^T. C ^A 


1.5 


frameshift 


9.0 
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Table 2 - Neutral evolution robustness and mutation data. 

Each row is for a different generation, T . Entries of NA indicate that no measurement was made. The {mnt) 
and {niaa) are the average number of nucleotide mutations and nonsynonymous mutations, respectively. 
Numbers in parentheses are total counts over the total samples. Subscripts indicate the population type: 
U for imselected, P for polymorphic, and M for monomorphic. For the unselected and monomorphic 
populations, numbers represent averages of all replicates. For the polymorphic population, numbers are for 
a random sample of functional mutants. {T)p and {T)m are the fraction of functional mutants out of 435 
assayed. 



T 


{m„t)u 


{maa)u 


{mnt)p 


{maa)p 































0.48 (210 / 435) 


0.48 (210 / 435) 


1 


NA 


NA 


NA 


NA 


0.1 (3 / 22) 


0.3 (6 / 22) 


0.48 (208 / 435) 


NA 


2 


NA 


NA 


NA 


NA 


0.4 (9 / 22) 


0.8 (17 / 22) 


0.49 (215 / 435) 


NA 


3 


5.0 (20 / 4) 


3.5 (14 / 4) 


2.7 (27 / 10) 


1.4 (14 / 10) 


1.0 (23 / 22) 


0.4 (9 / 22) 


0.49 (215 / 435) 


NA 


4 


NA 


NA 


NA 


NA 


1.5 (32 / 22) 


0.7 (15 / 22) 


0.48 (208 / 435) 


NA 


5 


NA 


NA 


NA 


NA 


2.2 (48 / 22) 


1.1 (25 / 22) 


0.45 (197 / 435) 


0.43 (185 / 435) 


6 


9.8 (39 / 4) 


7.5 (30 / 4) 


5.5 (55 / 10) 


2.1 (21 / 10) 


2.6 (58 / 22) 


1.4 (31 / 22) 


0.46 (198 / 435) 


NA 


7 


NA 


NA 


NA 


NA 


3.1 (69 / 22) 


1.8 (39 / 22) 


0.52 (227 / 435) 


NA 


8 


NA 


NA 


NA 


NA 


3.4 (74 / 22) 


1.8 (40 / 22) 


0.46 (200 / 435) 


NA 


9 


13.0 (52 / 4) 


10.3 (41 / 4) 


6.7 (61 / 9) 


3.1 (28 / 9) 


3.7 (82 / 22) 


2.1 (46 / 22) 


0.47 (203 / 435) 


NA 


10 


NA 


N.'\ 


N.A. 


NA 


4.2 (92 / 22) 


2.4 (52 / 22) 


0.46 (199 / 435) 


0.40 (175 / 435) 


11 


NA 


N.V 


X,\ 


N.V 


4,(i (11)2 / 22) 


2,5 (.">(! ,' 22) 


11,48 (21)7 4:!.".) 


X,\ 


12 


llj.8 (Ij7 / 4J 


12,.''. {TAI I 4) 


7,,'i (7l( ,; <J) 


(-'M : <J) 


4,9 (iU7 / 22) 


2,(1 158 / 22J 


0,52 (228 /■ 4,35 ) 


NA 


13 


NA 


NA 


NA 


NA 


5.0 (110 / 22) 


2.7 (60 / 22) 


0.52 (227 / 435) 


N.A. 


14 


NA 


NA 


NA 


NA 


5.3 (116 / 22) 


2.9 (64 / 22) 


0.50 (216 / 435) 


NA 


15 


NA 


NA 


10.3 (227 / 22) 


3.8 (83 / 22) 


5.6 (123 / 22) 


3.0 (67 / 22) 


0.50 (219 / 435) 


0.39 (171 / 435) 


16 


NA 


NA 


NA 


NA 


5.8 (127 / 22) 


3.0 (67 / 22) 


NA 


NA 


17 


NA 


NA 


NA 


NA 


6.0 (133 / 22) 


3.1 (69 / 22) 


NA 


NA 


18 


NA 


NA 


NA 


NA 


6.3 (137 / 22) 


3.2 (71 / 22) 


NA 


NA 


19 


NA 


NA 


NA 


NA 


6.3 (138 / 22) 


3.3 (72 / 22) 


NA 


N.A. 


20 


NA 


NA 


NA 


NA 


6.6 (145 / 22) 


3.4 (75 / 22) 


NA 


0.37 (160 / 435) 


21 


NA 


NA 


NA 


NA 


6.9 (152 / 22) 


3.6 (79 / 22) 


NA 


NA 


22 


NA 


NA 


NA 


NA 


7.1 (156 / 22) 


3.7 (81 / 22) 


NA 


NA 


23 


NA 


NA 


NA 


NA 


7.2 (158 / 22) 


3.7 (81 / 22) 


NA 


NA 


24 


NA 


NA 


NA 


NA 


7.3 (161 / 22) 


3.8 (83 / 22) 


NA 


NA 


25 


NA 


NA 


NA 


NA 


7.7 (169 / 22) 


4.0 (87 / 22) 


NA 


0.39 (169 / 435) 
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Additional material 

Additional file 1 - Sequence of the parent P450 used start neutral evolution. 

FASTA file with sequence of the Rl-11 P450 BME used as the neutral evolution parent. This sequence was 
isolated after the equilibration evolution. 

Additional file 2 - Information about sequences from neutral evolution experiments. 

The entries give the name of the mutant, the number of nonsynonymous and nucleotide mutations relative 
to the Rl-11 parent, the [urcajso value if measured, the T50 value if measured, the percent of the parental 
expression level if measured, and then a list of all of the mutations. Amino acid mutations are numbered in 
the standard P450 numbering scheme. The names of the mutants indicate their origin. Names beginning 
with "P-G3" arc randomly chosen functional mutants from generation 3 of the polymorphic population, 
etc. Names of the form "PI," "P2,", etc. are the 22 functional mutants that were randomly chosen from 
the final (generation 15) polymorphic population. Numbers P5 and P12 are missing because two of the 
original 24 randomly selected polymorphic population members were randomly chosen to be discarded after 
it was discovered that two of the 24 monomorphic replicates were contaminated. Names beginning with 
"Ul" indicate that sequences are from the first unselected replicate, etc. Names beginning "Ml" indicate 
sequences are from the first monomorphic replicate, etc. Replicates "M9" and "MIO" were discarded due to 
contamination during the experiment. For each replicate, we sequenced each new functional mutant. The 
last functional mutant after 25 generations represents the final sequence for that replicate, and is given an 
abbreviated name without the generation suffix. 

Additional file 3 - Thermostability measurements. 

Raw data from the T50 thermostability measurements. 

Additional file 4 - Urea stability measurements. 

Raw data from the [urea] 50 thermostability measurements. 

Additional file 5 - Correlation of thermal and urea stabilities. 

The T50 and [urea] 50 values are highly correlated. 

Additional file 6 - Sequence of initial P450 used to start equilibration evolution. 

FASTA file with sequence of the 21B3 P450 BM3 heme domain described in [23]. This P450 was used as 
the initial parent to start the equilibration evolution. 

Additional file 7 - Mutations accumulated during equilibration evolution. 

The file lists the mutations in the 46 P450 variants selected at the end of the equilibration evolution. Each line 
gives the name of the variant, with the prefix indicating whether it came from the Rl or R2 population. The 
next entries give the number of nucleotide and nonsynonymous mutations. All of the individual mutations 
relative to 21B3 are then listed. Amino acid mutations are numbered in the standard P450 numbering 
scheme, with the threonine after the N-terminal methionine given the number one. 
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